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ABSTRACT 

We present Monte Carlo calculations of Lya radiative transfer through optically thick, spherically symmetric, 
collapsing gas clouds. These represent simplified models of proto-galaxies that are caught in the process of their 
assembly. Such galaxies produce Lya flux over an extended solid angle, either from a spatially extended Lya 
emissivity, or from scattering effects, or both. We present a detailed study of the effect of the gas distribution and 
kinematics, and of the Lya emissivity profile, on the emergent spectrum and surface brightness distribution. The 
emergent Lya spectrum is typically double-peaked and asymmetric. In practice, however, we find that energy 
transfer from the infalling gas to the Lya photons -together with a reduced escape probability for photons in the 
red wing- causes the blue peak to be significantly enhanced and the red peak, in most cases, to be undetectable. 
This results in an effective blueshift, which, combined with scattering in the intergalactic medium, will render 
extended Lya emission from collapsing protogalaxies difficult to detect beyond redshift z <; 4. We find that 
scattering flattens the surface brightness profile in clouds with large line center optical depths (tq > 10^). A strong 
wavelength dependence of the slope of the surface brightness distribution (with preferential flattening at the red 
side of the line) would be a robust indication that Lya photons are being generated (rather than just scattered) 
in a spatially extended, collapsing region around the galaxy. We also find that for self-ionized clouds whose 
effective Lya optical depth is < 10^, infall and outflow models can produce nearly identical spectra and surface 
brightness distributions, and are difficult to distinguish from one another The presence of deuterium with a cosmic 
abundance may produce a narrow but detectable dip in the spectra of systems with moderate hydrogen column 
densities, in the range 10'**- 10^" cm"^. Finally, we present a new analytic solution for the emerging Lya spectrum 
in the Umiting case of a static uniform sphere, extending previous solutions for static plane-parallel slabs. 

Subject headings: galaxies: formation - galaxies: halos - quasars: general - radiative transfer - cosmology: 
theory - intergalactic medium 



1. INTRODUCTION 

The early stages of galaxy formation are characterized by a 
spatially extended distribution of gas. As gas cools and settles 
inside its host dark matter halo, most of its gravitational binding 
energy is radiated away. In dark matter halos with Tvii > 10"^ 
K, the majority of this radiation can emerge in the Lya line 
(Haiman et al. 2000; Fardal et al. 2001). Protogalactic gas 
clouds may therefore be detected as spatially extended sources 
in the Lya line. In the presence of ionizing sources enclosed 
by the collapsing gas, the total Lya flux from these clouds may 
be significantly boosted, as photo-ionization and subsequent re- 
combination converts the ionizing continuum photons into Lya 
radiation. The ionizing sources may consist of young massive 
stars (which are hot enough to produce substantial ionizing ra- 
diation) and/or quasars. For example, Haiman & Rees (2001) 
have shown that in the presence of a luminous central quasar, 
this boost in the amount of Lya emission from a spatially ex- 
tended protogalactic gas cloud may be several orders of magni- 
tude. 

Based on the three Lya emission mechanisms sketched above 
- cooling radiation, or recombination powered by the ionizing 
radiation of stars or quasar black holes, - copious Lya emis- 
sion may be expected from the early stages of galaxy forma- 
tion, during which the gas is spatially extended and collapsing. 
Simple estimates of the total Lya line flux show that it can be 
detectable out to high redshifts with modern instruments. In- 
deed, Wide-field narrow-band surveys, such as for example 
the Large Area Lyman Alpha (LALA) survey (Rhoads et al. 
2000) have been very fruitful in finding young galaxies at red- 



shifts z = 4.5-6.5. Although the majority of the observed Lya 
emitters appear point-like (within the detection limits), several 
tens of extended Lya sources have currently been observed at 
z 3 (e.g Steidel et al. 2000; Bunker et al. 2003; Matsuda et 
al. 2004). One may expect that a gas cloud that is caught in 
the early stages of galaxy assembly to be relatively metal-free, 
since not much star-formation has taken place. The observa- 
tions of Weidinger et al. (2004, 2005) suggest that such a pris- 
tine (not polluted by metals) gas cloud around a z '--^ 3 quasar has 
already been found. A similar extended Lya "fuzz" has been 
discovered around 4 additional quasars recently by Christensen 
(2005). Barkana & Loeb (2003) have shown that evidence for 
gas collapse around quasars exists in the form of an absorption 
feature in the quasar's Lya emission line. 

Future searches for both Lya emitters and blobs will un- 
doubtedly reveal more Lya emitters, point-like and extended. 
Deeper observations and/or higher resolution spectra may con- 
tain detailed information on the gas distribution and kinematics. 
Therefore, it is timely to ask the following question: How can 
we use the Lya radiation itself, to identify these objects as col- 
lapsing protogalaxies-and to constrain the mechanisms power- 
ing their Lya emission ? 

Quantifying the properties of Lya radiation emerging from 
collapsing protogalaxies, such as its expected surface bright- 
ness profile and spectral energy distribution, is difficult because: 
(i) it is not well understood how and where the Lya photons are 
generated, (ii) the optical depth t in the line exceeds unity for 
columns of neutral hydrogen > 10'"^ cm so that protogalax- 
ies are expected to be optically thick with t 10^^, necessitating 
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computationally expensive radiative transfer calculations, and 
(iii) complex gas density and velocity distributions (e.g. Kunth 
et al. 1998), as well as the presence of dust (e.g. Neufeld 1991; 
Chariot & Fall 1993; Hansen & Oh 2006), complicate any sim- 
plified predictions. Additionally, at redshifts z > 3, the inter- 
galactic medium (IGM) is optically thick at the Lya frequency, 
and wiU modify the observed spectrum. 

In this paper, we attack this problem with Monte Carlo com- 
putations, similar to those used in previous papers (e.g. Ahn 
et al. 2000, 2002; Zheng & Miralda-Escude 2002; Cantalupo 
et al. 2005) of Lya transfer through a range of models that 
represent collapsing protogalactic gas clouds. The main goal 
of this work is to gain basic insights into the properties of the 
emerging Lya radiation, and into how these properties depend 
on the assumed gas density distribution, gas kinematics and 
on where the Lya is generated. Furthermore, we study how 
these predictions change in the presence of an embedded ioniz- 
ing source. Our focus is on the in the earliest stages of galaxy 
formation, when there is still significant gas infall with an ex- 
tended spatial distribution. Therefore the geometry we study 
is that of a spherically symmetric, collapsing gas cloud, which 
is obviously oversimplified. However, the simple geometry is 
computationally less expensive and allows us to explore a wide 
range of Lya emissivity profiles and gas distributions that may 
be appropriate for the collapsing protogalaxies. Additionally, 
the absence of complicated geometrical effects makes our re- 
sults easier to interpret and may fulfill the same role for future 
(and existing, e.g. Tasitsiomi 2005), more sophisticated 3-D 
simulations, that analytic solutions have played for the present 
work. Furthermore, we will argue that several of our result are 
unlikely to change for more realistic, complex gas distributions 
and kinematics. In a companion paper (Dijkstra et al. 2006, 
hereafter Paper II), we present a detailed comparison of our 
models with several existing observations of Lya emitters. 

The rest of this paper is organized as follows. In § 2, we 
briefly review the basics of Lya transfer. In § 3, we describe 
our code, along with several tests we performed to check its ac- 
curacy. In § 4, we describe how deuterium was included in our 
calculations, and a motivation for why and when this is relevant. 
In § 5, we apply our code to Lya transfer calculations of neutral 
collapsing gas clouds (without any embedded ionizing sources). 
In § 6, we present present our main results on the emergent 
spectra and surface brightness profiles. In § 7, we investigate 
the additional effect of an embedded ionizing source. In § 8, we 
discuss the impact of the IGM and various caveats in our calcu- 
lations. Furthermore, we compare our results with other work. 
Finally, in § 9, we present our conclusions and summarize the 
impUcations of this work. The parameters for the background 
cosmology used throughout this paper are $7,,, = 0.3, = 0.7, 
Qb = 0.044, h = 0.7, based on Spergel et al. (2003). 

2. LYa RADIATIVE TRANSFER BASICS 

The basics of transfer of Lya resonance line radiation has 
been the subject of research for many decades and is well un- 
derstood (e.g. Zanstra 1949; Unno 1952; Field 1959; Adams 
1972; Harrington 1973; Neufeld 1990). The goal of this sec- 
tion is to review the basics of Lya scattering that is relevant to 
understanding this paper. It is convenient to express frequencies 
v in terms of x = (u-i/q)/ Auo, where Auo = Vihfo/c, and Vih 
is the thermal velocity of the hydrogen atoms in the gas, given 
by Vth = ^yiksT /nip, where kg is the Boltzmann constant, T 
the gas temperature, nip the proton mass and uq is the central 



Lya frequency, i^o = 2.47 x 10 Hz. 

The optical depth through a column of hydrogen, Nm, for a 
photon in the line center, tq is given by 

To =Nm(To=Nmfi2-;;^^^ = 

= 8-3xlO''(,^^,i^)(^)""', (1) 

where ctq is the Lya absorption cross section in the line center, 
/i2 = 0.4167, the oscillator strength and nie and e are the mass 
and charge of the electron, respectively. The optical depth for a 
photon at frequency x is reduced to 

I^=H(a,x)=^ r-^^, = l (2) 
To (y-x)^+a^ I ^mg, 

where a is the Voigt parameter and is the ratio of the Doppler 
to natural line width, given by a = A2i/4'kAi^d = 4.7 x 10"^ 
(13 km /vti,), where A21 is the Einstein A-coefficient for the 
transition. The transition between 'wing' and 'core' occurs at 
X ~ 2.5 -4.0 for a ~ lO'^-lO'^ 

When a Lya photon of frequency x, is absorbed by an atom, 
it re-emits a Lya photon of the same frequency in its own 
frame, unless the atom is perturbed while in the excited 2p 
state. These perturbations include collisions with an electron 
and photo-excitation or ionization by another photon. In this 
case, the energy of the re-emitted Lya photon, Xo, is not equal 
to X,. In most astrophysical conditions these perturbations can 
be neglected (see appendix D), rendering the scattering coher- 
ent in the frame of the atom. 

Due to the atom's motion however, the scattering is not co- 
herent in the observer's frame. A simple relation exists between 
the incoming and outgoing frequency, x, and Xo, and the veloc- 
ity of the atom, Va- 

Xo=Xi -I- +g(krko-l) + 0{\ — , (3) 

Vth Vth V- c / 

where ki and k,, are the propagation directions of the incom- 
ing and outgoing photons, respectively. The second to last term 
on the RHS is the 'recoil' term, which accounts for the trans- 
fer of a small amount of momentum from the photon to the 
atom during each scattering (exactly as in Compton scattering, 
e.g. Rybicki & Lightman 1979, p. 196). The factor g is the av- 
erage fractional amount of energy transferred per scattering. 
Field (1959) showed that g can be written as ^ = hAvo/lkT 

= 2.6 X lO"'' (13 km s~'/vth). The effect of recoil is (and has 
been verified to be) negligible in the appUcations presented in 
this paper (Adams 1971). 

Because x, and x,, are related via eq. (3) this case is referred to 
as 'partially coherent'. The redistribution function, q(xo,Xi)dxo, 
gives the probability that a photon that was absorbed at fre- 
quency X, is re-emitted in the frequency range Xg ± dx„/2 (in 
the observer's frame). Under the assumption that the velocity 
distribution of atoms is (locally) Maxwellian, Unno (1952) and 
Hummer (1962) calculated the redistribution function analyti- 
cally for partially coherent scattering ( their ^//(xo,x,)). Their 
calculation provides a good test case for our code, which is dis- 
cussed in Appendix B. 

Photons in the core have a short mean free path, and the 
rate at which they are scattered is high. The majority of these 
scattering events leave the photon in the core, but a chance en- 
counter with a rapidly moving atom may move it to the wing. 
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where its mean free path is much larger. In moderately optically 
thick media, tq ^ 10^, one such encounter may be enough to en- 
able the photon to escape in a single flight. According to eq. (1) 
this optical depth corresponds to a column density of 2.4 x lO'^ 
cm"^ for a gas at T = 2 x 10"* K. This is approximately the col- 
umn density of a Lyman limit system. Adams (1972) demon- 
strated that in the case of extremely opticaUy thick media, tq > 
10^ /a, photons do not escape in a single flight, but rather in 
a single 'excursion' . During this excursion the photon scatters 
multiple times in the wing, and diffuses through real and fre- 
quency space. Because of the significantly increased mean free 
path, a photon predominantly diffuses spatially, while it is in 
the wing. Each scattering in the wing pushes the photon back 
towards the core by an average amount of -1/x (Osterbrock 
1962). 

The transfer of Lya photons through a static, uniform, non- 
absorbing, plane-parallel scattering medium (hereafter referred 
to as "slab") has been studied extensively. The spectrum of Lya 
photons emerging form both moderately and extremely opti- 
cally thick slabs is double peaked and, when bulk motions in the 
gas are absent and when recoil is ignored, symmetric around the 
origin. The spectra of Lya photons emerging from extremely 
optically thick slabs have been calculated analytically. For the 
case in which photons are injected in the center of the slab, 
in the Une center, the location of the peaks at ±Xp is given by 
Xp = 1.06(aro)'/^ (Harrington 1973; Neufeld 1990), in which tq 
is the total line center optical depth from the center to the edge 
of the slab ^ 

3. THE CODE AND TESTS 

3.1. The Code 

Our Monte Carlo computations are based on the method pre- 
sented by ZM02. Our code focuses on spherically symmetric 
gas configurations. We use A^grid spherical shells, each with their 
own radius, density, neutral fraction of hydrogen and velocity. 
The temperature is assumed to be the same in each shell, al- 
though the code can easily be adjusted to account for tempera- 
ture variations. Throughout this paper, we used A'grid = 1000. 
Tests performed with 10 times more shells produced identi- 
cal results. Runs performed with ^ 10 times less shells still 
produced visually indistinguishable results. Because runs with 
A'grid = 100 barely reduced the computing time, A^grid = 1000 was 
used. The shells were spaced such that every shell contains 
roughly an equal column of hydrogen from its inner to its outer 
radius. 

The creation and transfer of individual Lya photons is de- 
scribed below. In the transfer problem, we identify three rel- 
evant frames: the 'central observer' frame; the frame of the 
shell through which the photon is traveling; and the frame of the 
atom that scatters the Lya photon (the frame of the atom differs 
from that of the shell it is in, because of the atom's thermal mo- 
tion). The 'central observer' frame is the frame of an observer 
located in the center of the sphere. The reason that the observer 
is put in the center is that it fully exploits spherical symmetry. 
The bulk motion of the gas is simply pointed toward (for infall) 
or away (for outflows) from the observer. Photons propagat- 
ing outwards and inwards, propagate away and towards the ob- 
server, respectively. Form now on the 'central observer' frame 
will be referred to simply as the observer's frame. 

'Note that the value of Xp quoted by Harrington (1973) and Neufeld (1990), 
is actually Xp = 0.88(aro)'^^. This is due to their different definition of H{a,x), 
which is lower by a factor of (Ahn et al. 2002). 



A quantity '2' measured in the frame of a shell, observer 
and the atom is denoted by Q', Q and Q, respectively. /?„, with 
n = 1 , 2 , 3 , . . ., are random numbers generated between and 1 . 

Step 1. Given the emissivity per unit volume as a function 
of radius, j{r), we generate the radius of emission r^m from the 
probabiUty distribution function: 



-Jir')dr' 



(4) 



where rmax is the radius of the gas sphere's edge. The loca- 
tion of the photon in the sphere, r = ipx,Py,Pz)7 is given by 
{Pi:Py,Pz) = rem (S2 C3, ^2 ^3, C2), where C„ and 5„, with n = 
1 , 2 , 3 , . . . , are cos 27r/?„ and sin 2nR„ , respectively. The emissiv- 
ity jir) depends on what model we consider, which is discussed 
in § 5.2. 

Step 2. Given the location of emission, we generate the di- 
rection of emission ki = {kx,ky,k^={R4 C5, Ss, (1 -/?|)°^). In 
the case the gas has a bulk velocity, the frequency of the emit- 
ted photon is given by x, = x'j + kj ■ v'buik/i'th, where Vbuik is the 
bulk velocity vector of the gas at the location of emission. For 
simplicity, we set x[ = 0, which means that the Lya is initially 
always injected at the line center in the frame of the gas sheU in 
which it is emitted. In reality x\ is distributed according to the 
Voigt profile (eq. 2). This minor simplification does not affect 
our results at all. The reason is that the Lya photon is redis- 
tributed in frequency each time it is scattered, rapidly erasing 
any 'memory' of its initial frequency. For our spherically sym- 
metric models, Vbuik = Vbuik'71'1' where Vbuik is the infall/outflow 
speed of the gas at radius r. The sign of Vbuik is negative/positive 
for infall/outflow, respectively. 

Step 3. The optical depth r the photon is allowed to travel is 
given by -ln[/?6]- To convert r to a physical distance, we calcu- 
late the physical distance to the edge of the shell in which the 
photon is emitted in the direction ki. The corresponding optical 
depth to the edge of the shell is = de «h <7() 4>i?^i), where hh 
is the number density of neutral hydrogen atoms in the shell. 
If Te > T, then the new location is given by F+ Xki, where 
\ = de r/Te- li Te < T then, we obtain the frequency of the 
photon in the frame of the next sheU from x- = x, - ki ■ i^buik/vth, 
where Vbuik is now the bulk velocity vector of the gas where the 
photon penetrates the next shell. In the next shell r is replaced 
by T- r<.. This process is repeated until the photon has traveled 
the generated optical depth r to its new position f. 

Step 4. When |r| > r„ax , the photon's frequency and the 
angle under which it passed through the surface of the sphere 
are recorded. Otherwise, the photon is scattered by a hydrogen 
atom, in which case we proceed to the next step. 

Step 5. Once the location of scattering, r, has been deter- 
mined, the total velocity of the atom, Va, that scatters the Lya 
photon is generated. The total velocity is given by the sum 
of the bulk velocity of the gas plus the thermal velocity: Va = 
Vbuik + The thermal velocity, Vth, is decomposed into two 
parts. Following ZM02, the magnitude of the atom's velocity 
in the direction ki, M||, is generated from eq. (5), 

rR9 „-/ 



M|| = 



,dy. 



(5) 



TTH(a,4)J_^ {x[-yf+a^ 
This distribution reflects the strong preference for photons at 
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frequency x-, to be scattered by atoms to which they appear ex- 
actly at resonance, uu = x'^ . For large x-, however, the number of 
these atoms reduces as cx e and the photon is Ukely scattered 
by an atom to which it appears far in the wing. 

The two (mutually orthogonal) components perpendicular to 
ki, MX 12, are drawn from a Gaussian distribution with zero mean 
and standard deviation 2"^/^ v,/,. Following the procedure de- 
scribed in Numerical Recipes (Press et al. 1992) 

M_Li = ^-\n[Rn] cos(27r/?io) 

M±2 = \/-\n[R\i\ sin(27r/?io), (6) 

where m_li,2 is in units of Vth- At this particular step in the code, 
it is possible to speed up the code significantly, which is de- 
scribed below (§ 3.1.1). 

Step 6. In the frame of the atom, the distribution of the direc- 
tion of the re-emitted photon with respect to that of the incom- 
ing photon is given by a dipole distribution. 



Ri = - 



(I + ^'Va^', 



(7) 



where ii = ki-ko. Given ki, we have = iiki+ ( 1 - iJ,^)k±, where 
k± is an arbitrary unit vector perpendicular to ki. kg is calcu- 
lated from ko via a Lorentz transformation into the observer's 
frame. Once Va and kg are known, can be calculated using 
eq. 3. After replacing x,- and ki by Xg and kg, respectively, we 
return to step 3. 

3.1.1. Accelerated Scheme 

It is possible to speed up the code using a method described 
by Ahn et al. (2002). A critical frequency, Xcrit, defines a tran- 
sition from 'core' to 'wing'. Once a photon is in the core 
(|x'| < Xcrit), it is possible to force it (back) into the wing, by 
only allowing the photon to scatter off a rapidly moving atom. 
This is achieved by drawing mxi and u±^2 {Step 5. in § 3.1) 
from a truncated Gaussian distribution, which is Gaussian for 
u > Xcrit, but for M < Xcrit- To be more precise: 



M_L1 = 



M±2 ■■ 



- . 



-hi[/?ii] cos(27r/?io) 



-ln[/?ii] sin(27r/?io), 



(8) 



where u_li,2 is again in units of Vth. 

The reason it is allowed to skip core scattering is that spatial 
diffusion occurs predominantly in the wings (see § 2). Skip- 
ping scattering events in the core effectively sets the mean free 
path of core photons to 0. Clearly, Xcrit = restores the non- 
accelerated version of the code. Unless stated otherwise, in 
this paper we used Xcrit = 3. We found that larger values of Xcnt 
changed the appearance of the spectrum, especially the pres- 
ence of the red peak (see below) was strongly affected. 

3.2. The Tests 

We performed various tests on the code, three of which are 
described in Appendix B. The other test is explicitly described 
here, since it involves a new analytic solution we derived in 
Appendix C. The test is similar to one generally done by other 
groups (e.g. Ahn et al. 2000, 2002; Cantalupo et al. 2005, ZM02), 
but differs in the geometry of the gas distribution. Lya photons 
are injected at x = in the center of a uniform, non-absorbing. 



static sphere with the gas at T = 10 , which corresponds to 
a= 1.5 X 10"^. The line center optical depth from the center 
to the edge of the sphere is tq. In the upper left panel of Fig- 
ure 1 the emergent spectra are shown for three different values 
of T(): T() = 10^, lO*" and 10^ and compared to the analytic solu- 
tion given below. 



J{x) = 



24aTo 



, 1+cosh 



27 aTo 



(9) 



The derivation of this analytic solution for a sphere is com- 
pletely analogous to the derivation given by Harrington (1973) 
and Neufeld (1990), and can be found in Appendix C. As men- 
tioned in § 2, the emergent spectra are double-peaked and sym- 
metric around the origin, with Xp increasing with tq. For a given 
value of To, Xp = 0.92(aro)'/^ (vs. Xp = 1.06(aTo)'''^ for a slab). 
The spectra are normaUzed (as is the analytic solution) such that 
the area under all curves is (27r)~' . The figure shows that the 
agreement is good. In Appendix B other tests are performed on 
the code that involve core scattering and a case in which the gas 
has bulk motions. Here, we briefly sunmiarize these tests: In 
the upper right panel the redistribution functions, qnix',x), for 
X = (blue), X = 2 (red) and x = 5 (green), as extracted from the 
code are compared with the analytic solutions given in Hum- 
mer (1962) and Lee (1974) (solid lines). The total number of 
scatterings a Lya photon undergoes before it escapes from a 
slab of optical thickness 2to as extracted from the code (cir- 
cles) is shown in the lower left panel. Overplotted as the red- 
solid line is the theoretical prediction by Harrington (1973). In 
the lower right panel we show the spectrum emerging from an 
infinitely large object that undergoes Hubble expansion. The 
histogram is the output from our code, while the green-solid 
line is the (shghtly modified) solution obtained by Loeb & Ry- 
bicki (1999, hereafter LR99) using their Monte Carlo code. As 
Figure 1 shows, the code passes all these tests well. 

4. DEUTERIUM 

Although it has a low abundance, accounting for deuterium 
may be relevant when predicting spectra emerging from Lya 
emitters. The resonance frequency of the deuterium hya ab- 
sorption line is blueshifted by 82 km s"^ with respect to that 
of hydrogen. For a gas at T = 74 x 10^ K, this corresponds to 
xd = +6.3 (74)"^/^. According to eq. (2), the hydrogen optical 

depth at this frequency is reduced to Txo,h ~ 6.6 x 10~^ tI^^tq. 
At this frequency, the photon is (by definition) exactly in the 
deuterium line center, and the optical depth due to deuterium is 
rxD,D = 2^I\D/H]to ~ 4.4 x lO'Vo, where [D/H] is the num- 
ber of deuterium atoms per hydrogen atom^, [D/H] = 3 x 10"^ 
(Buries & Tytler 1998). The contribution of deuterium to the 
optical depth, / tq, is shown as a function of frequency in 
the range x = 3 - 10 (red-solid line in Figure 2). The black- 
dotted line) shows the original Txa/tq given by eq. (2). The 
(blue-dashed line) shows the sum of the two. The gas is as- 
sumed to be at °K . Clearly there is a range in frequencies, 
X ~ 5 - 8, where deuterium dominates the optical depth and reg- 
ulates the transfer of Lya photons. Despite extensive research 

^The factor 2'^^ reflects that deuterium is two times heavier than hydrogen, 
which reduces its thermal velocity by 2^/^, which according to eq. (1) increases 
the line center optical depth by 2'/^. Additionally, the reduced thermal velocity 
of deuterium raises its Voigt parameter by 2'/^. 
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Fig. 1 . — The panels in this figure show four different tests of our Monte Carlo code against known analytic solutions and Monte Carlo calculations by Loeb & 
Rybicki (1999). Upper left panel: The spectrum emerging from a uniform spherical gas cloud, in which Lya photons are injected in the line center, .v = 0. The total 
line center optical depth, tq from the center to the edge is 10^ and 10^ for the narrowest {blue) and broadest (green) profile, respectively. The intermediate spectrum 
(red) corresponds to tq = 10^. Overplotted as the black dotted line is the theoretically derived spectrum, eq. (9). The tests that are shown in the other panels are 
described in Appendix B. Upper right panel: The colored histograms are the redistribution functions, qii(x' ,x), for jc = (blue), x = 2 (red) and ,r = 5 (green), as 
extracted from the code. The sohd lines are the analytic solutions given in Hummer (1962) and Lee (1974). Lower left panel: The total number of scatterings a 
LyQ photon undergoes before it escapes from a slab of optical thickness 2tq as extracted from the code (circles). Overplotted as the red— solid line is the theoretical 
prediction by Harrington (1973). Lower right panel:. The spectrum emerging from an infinitely large object that undergoes Hubble expansion. The histogram is the 
output from our code, while the green-solid line is the (slightly modified) solution obtained by Loeb & Rybicki (1999) using their Monte Carlo algorithm. 



on Lya radiative transfer, this detail has been neglected, to our 
knowledge, in all previous studies. 

To include deuterium in our code we added the contribu- 
tion of deuterium, Tj £>, to Tj. In addition to this, a step was 
added to the code following Step 4. (§3.1), that determines 
whether a scattering occurs by a hydrogen or deuterium atom. 
A random number Ro was generated between and 1. Let 
Ph = Tx.h/{tx,h + Tx.d\ be the probability that the Lya photon 
is scattered by hydrogen. When Ro < Ph, the scattering is done 
by hydrogen, in which we proceed as described under Step 5. 
in § 3.1. Otherwise, the scattering is done by deuterium and its 
velocity vector is generated in exactly the same way, but with 
a and x replaced by 2'/^a and 2'''^[x-82 km s"'/v,/,], respec- 
tively. 

Figure 3 demonstrates the impact of deuterium on the Lya 
spectrum emerging from a uniform sphere in which the pho- 
tons are injected at line center (x = 0) in the sphere's center 
The gas is at T = lO'* °K and the total column of hydrogen 
from the center to the edge of the sphere is Nhi = 1.2 x 10'^ 



cm"^, yielding tq = 7.3 x 10^ (eq. 1). The emerging spectrum is 
shown with ired-soUd-histogram) and without iblack-dotted- 
histogram) deuterium taken into account. The blue-dashed line 
is the theoretically predicted spectrum, eq (9). The reason the 
theoretical curve does not fit very well is that the product aro = 
343, while Neufeld (1990) requires otq to be aro > lO"' for the 
theoretical curve to provide a good fit. The imprint deuterium 
leaves on the spectrum is strong. The main reason is that the 
presence of deuterium reduces the mean free path for photons 
near x = xq, which reduces the probability that they escape the 
medium at this frequency. We find that even for a 10 times 
lower deuterium abundance, the imprint in the spectrum is still 
noticeable. 

For deuterium to leave an observable imprint in the spec- 
trum, it should affect the frequencies at which a non-negligible 
fraction of Lya photons emerges (as in Fig. 3). For 7^=1, for 
example, we expect deuterium to leave a spectral feature for 
columns in the range A^hi = 10'^- 10^° cm"^, since this range 
would translate to Xp ~ (aro)'''"' ^3-12. It is interesting that 
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Fig. 2. — The contribution of a cosmological abundance of deuterium 
{[D/H] = 3 X 10"^) to the Lya opacity for a gas at T = lO'* K. The ratio of the 
optical depth at frequency x and frequency .v = 0, Tj/^/tq = H{a,x), through 
static gas is shown. The black-dotted line is for hydrogen only ( and is given 
by eq. 2); the red-solid line is for deuterium only; the blue— dashed line is 
their sum. The figure shows that in the range x = 5 — 8, deuterium dominates 
the opacity. For lower/higher temperatures, this spectral modification occurs 
at higher/lower values of x with a relatively larger/smaller contribution from 
deuterium, respectively. 
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Fig. 3. — The effect of a cosmological abundance of deuterium {[D/H] - 3x 
10"^) on the spectrum emerging from a sphere with a center-to-edge column 
density of hydrogen oi Nm = 1.2 X lO" cm"^. The reduced mean free path of 
Ly a photons in the deuterium line center inhibits their escape, which modifies 
the spectrum. For this plot, the gas temperature was assumed to be 10'* K. 



this corresponds to the range of column densities in Lyman limit 
systems. Because the columns studied in this paper are much 
larger, typically 10^' - 10^^ cm"^, the spectral feature due to 
deuterium is much less prominent. Furthermore, bulk motions 
in the gas will most likely smear the absorption feature over a 
wide frequency range, which makes it even harder to observe in 



practise. A further challenge is that the frequency range of the 
deuterium feature corresponds to 3 thermal widths, or v ~ 40 
km s"' . To resolve this narrow feature, spectra are required with 
R ~ 20000. This resolution can be achieved with the High Res- 
olution Echelle Spectrometer (HIRES) on Keck. Its detectabil- 
ity depends, however, on the total flux of Lya emission. The 



brightest known Lya blobs are ~ 10 ergs s" 



-I -2 

cm 



(Steidel 



et al. 2000), which is 2 orders of magnitude fainter than the 
sources successfully observed with this instrument. We plan to 
address the feasibility of detecting deuterium in more detail in 
future work. 

5. LYa EMISSION FROM NEUTRAL COLLAPSING GAS 
CLOUDS WITHOUT STARS OR QUASARS 

5.1. Physics of Cloud Collapse and Lya Generation 

For gas to collapse inside its host dark matter halo, it must 
lose its gravitational binding energy. In the standard picture 
this binding energy is converted to kinetic energy, which is con- 
verted to thermal energy when the gas reaches the virial radius, 
where it is shocked to the virial temperature of the dark mat- 
ter halo, Tvir. When the cooling time is less than the dynam- 
ical time, the gas practically cools down to 10'* °K instantly. 
Since this gas is no longer pressure supported, it falls inward 
on a dynamical timescale (Rees & Ostriker 1977). Haiman et 
al. (2000) argued that a large fraction of the cooling radiation 
is emitted in the Lya line (see also the early work by Katz & 
Gunn 1991, who computed the Lya cooling radiation in a hy- 
drodynamical simulation). Birnboim & Dekel (2003) presented 
a refined version of the previous scenario, and used analytical 
arguments in combination with a 1-D hydro simulation to show 
that virial shocks do not exist in halos with M <^ 2 x IO^Mq 
(which at z = 3 translates to Vdrc ^130 km s"'). Consequently, 
the gas in these halos never gets shocked to the virial tempera- 
ture and the gas is accreted 'cold' (also see Thoul & Weinberg 
1995). 

The details of cloud collapse and the associated Lya emis- 
sion are unclear In Birnboim & Dekel (2003)'s scenario, the 
gas falls unobstructed into the center, where its kinetic energy 
is converted to heat, which is radiated away in soft X-Rays. As 
a result, an Hll region forms in which recombination converts 
~ 68% of the ionizing radiation into Lya. The details of the gas 
distribution in the center and on the mean energy of the ionizing 
photons determine the radius of this Hll region, and thus how 
centrally concentrated the Lya emitting region is. 

Haiman et al. (2000) describe an alternative scenario that is 
supported by 3-D simulations, that show that any halo is build 
up of mergers of pre-existing clumps. The denser colder clumps 
fall inwards, while experiencing weak secondary shocks from 
supersonic encounters with other density inhomogeneities. As 
a result, the clumps are continuously heated, which is balanced 
by cooling through Lya emission. In the case the clumps fall 
in at a constant speed, all the gravitational binding energy of 
a clump is converted into heat which is radiated away in the 
Lya line. The Lya emissivity associated with this scenario is 
spatially extended. 

The main difference in the above models is the degree of con- 
centration of the Lya emissivity. While we cannot realistically 
model the emissivity profiles ab-initio, we will below consider 
a wide range of profiles that should bracket the possibilities dis- 
cussed above. 

5.2. Modeling the Gas Kinematics and Lya Emissivity Profile 
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Below, we define our model explicitly, and list the model pa- 
rameters. We discuss the uncertainties of each parameter, and 
motivate the range of values we consider in each case. Input 
that is kept fixed in all models is the density distribution of the 
dark matter, poi/), which is given by a NFW profile with con- 
centration parameter c = 5. The gas temperature is assumed to 
be 10"' K for all models. As discussed in § 8.2, changing the gas 
temperature does not affect our main results. The density distri- 
bution of the baryons, pb{r), is assumed to trace the dark matter 
at large radii, with a thermal core at r < iRs/A (see Mailer & 
Bullock 2004, their eq. 9-11, where Rs =/?vir/c). Although 
the results are not presented in this paper, we found the effect 
of varying the concentration parameter c to have no significant 
impact on our results. 

The total (dark matter + baryons) mass range of the collaps- 
ing cloud, Mtot, covered in this paper is lO''^- IO'^Mq. A re- 
lated quantity is the circular velocity, Vdrc- The upper limit on 
this mass range is set by simple cooling arguments: the cool- 
ing time of objects exceeding lO'^M© exceeds the Hubble time 
(e.g Blumenthal et al. 1984). The lower limit corresponds to 
the mass that can collapse even in the presence of an ionizing 
background (lO'^M© corresponds to Vcirc ~ 50 km s~^ at z = 3, 
see e.g. Thoul & Weinberg 1996; Kitayama & Ikeuchi 2000; 
Dijkstra et al. 2004). The redshift range at which the system 
virializes, Zvk covered in this paper is Zvir = 3-6. The lower 
Umit is set by the redshift of the observed Lya blobs which are 
candidates for coohng radiation. The upper limit is set by the 
opacity of the IGM which drastically increases beyond z > 6 
(e.g. Fan et al. 2002). Other input that is varied includes: 

The velocity field, Vbuik('"), is given by a simple power law 
Vbuik('') = Vamp[''/''yi>]"- The vclocity profile is uncertain; it can 
conceivably increase or decrease with radius. A spherical top- 
hat model would have a = 1 and Vamp = Vciic-^ In more reaUs- 
tic initial density profiles the mean density within radius r de- 
creases smoothly with radius, causing the inner shells to be de- 
celerated more relative to the overall expansion of the universe. 
This enhances the infall speed at small r relative to that in the 
top-hat model. The result of this is a flatter velocity profile, 
i.e a < 1 . Accretion of massless shells onto a point mass re- 
sults in a = -1 /2 (Bertschinger 1985). Since it spans the range 
we can motivate physically, we study values of a in the range 
a e [-0.5, 1], but note a may in fact be ^ 0.5. For cases with 
a < 0, the velocity diverges as r ^ 0. This artificially boosts the 
emission at small r. To prevent this, these cases are modified to 
flatten at r < Tv: 



Vbulk('-) = V; 



amp 



l\. + Tvir 



(10) 



When r > r^, Vbuik = Vamp['"/ '"virl" and when r ^ 0, then Vbuik ^ 
Vmax- We set Vmax = 2vamp, wWch uniquely determines r^. The 
amplitude of the velocity field Vamp is also varied between and 

2 Vcirc* 

The emissivity, jir), which we take to be either 'extended' 
or 'central. The extended case refers to the scenario sketched 
by Haiman et al. (2000), in which the gas continuously emits a 
fraction of the change of its gravitational binding energy in 
the Lya line. The emissivity as a function of radius jir), in ergs 
s~' cm"-', is then given by 

^This can be seen by time-reversing solutions for trajectories initially follow- 
ing Hubble expansion, v oc r. 



G Mtot(< r) pb dr 
r2 dt' 



(11) 



where Mtot(< r) is the total mass enclosed within radius r, pt 
the mass density of baryons at radius r and dr/dt is simply Vhuik, 
the speed at which the gas at radius r is falling in. We assume 
/q,=1, which reflects extreme cases in which aU gravitational 
binding energy is converted into Lya. The results given below 
scale linearly with /q,. 

The 'central', jir) = lC6(r), models represent the extreme ver- 
sion of the scenario sketched by Birnboim & Dekel (2003), in 
which all Lya is only generated in the center. Alternatively, 
these central models may represent any model in which a cen- 
tral Lya emitting source is surrounded by a neutral collapsing 
gas cloud. We determine /C from we the total Lya luminosity, 
given by Liy^ = J dV j{r). For simpUcity, we set Liya equal 
in the extended and central models with otherwise identical pa- 
rameters. 

A brief summary of all models studied in this paper is given 
in Table 2 in Appendix A. All results below are obtained by 
setting Tniax = ''vir (§3.1), which imphes that we ignore the Lya 
transfer beyond the virial radius. The impact of the IGM sur- 
rounding the virialized halo can vary from location to location 
(e.g. due to stochastic density and ionizing background fluctu- 
ations and also from the proximity effect near bright QSOs). 
We take the approach of studying the effects of scattering in 
the viriahzed halo and in the IGM separately. In this paper, we 
only compute the Lya spectra as they emerge at the surface of 
the halos at ryir with a brief discussion of the effect of the IGM 
in § 8.1. the IGM is discussed in more detail in Paperll. 

6. RESULTS 

Before discussing the spectral features and surface bright- 
ness profiles, it is useful to express the total Lya luminosity, 
LLya ihya = J dV j{r), § 5.2), as a function of our model pa- 
rameters: 



^Lyo 
T042 



" A^^ ( !^lot\-^ (Vamp\ (\+Z^ir\-- (2- a\^'^ 



ergs 
sec ' 



(12) 

The scalings of L\ya with Mtot and Vamp and (1 +Zvir) are exact. 
The a dependence is a fit, that is accurate to within 10%, over 
the range a e [-0.5, 1]. Note that equation (12) can be written 
as: 



-^Lya 



= 0.66 



ire \'' ( ^^fnp \ ( ^ 
m S"V V Vcirc / V 1.' 



1.2 



ergs 
sec 



(13) 



10^2 ' V 116 km s-iy V Vcirc >' V 1.75. 
The dependence of LLya on ^drc can be understood from the 
following scaling relations: LLya ocU /t (U is the total gravita- 
tional binding energy and f is the timescale over which it is re- 
leased, which is taken to be the dynamical timescale of the halo) 

oc f//(/?vir/Vciic) = Vcirct//^vir OC VdrcM^/Ryir OC V^^^, where WC 

used fliat U oc GM^/R^u and v^^^ = GM/R^i,. 

In § 6.1 we show some examples of Lya spectra and sur- 
face brightness profiles emerging from several models. The 
main purpose of this is to gain basic insight into the Lya ra- 
diation emerging from neutral, collapsing gas clouds with no 
star formation or quasar activity. In § 6.2 the dependence of 
these properties on the model parameters is studied in more de- 
tail. Our main findings are summarized and and qualitatively 
explained in § 6.3. 
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Fig. 4. — Properties of tlie hya emission emerging from our four fiducial models (models 1., 2, 1 c and 2 c) of neutral collapsing gas clouds. The left panels 
show spectra in which the lower horizontal axis shows the normalized frequency variable .v. The upper horizontal axis shows what velocity offset this frequency 
corresponds to. Values of .v > correspond to higher frequencies, or lower redshifts, which corresponds to lower recessional velocities. The left vertical axis shows 
the normalized flux, while the right vertical axis has the flux density in /iJy. Processing by the IGM is ignored in the plot. The right panels show the corresponding 
surface brightness distribution. The thick dotted horizontal line corresponds to the detection limit in a recent survey for Lya blobs by Matsuda et al. (2004). The 
model number is shown in the top right comer. The details of each model can be found in Table 2. The black-dashed line, in the upper P,vo panels in the right 
column, is the surface brightness profile that would emerge, in the absence of scattering. For the central models, this line corresponds to a <5-function at 0", and is 
therefore not shown. For all infall models, the Lya line emerges with a net blueshift. 



6.1. Examples of Intrinsic Spectra and Surface Brightness 

Profiles 

We show the emerging spectrum {left panels) and surface 
brightness profile {right panels) for models 1., 1 c, 2 and 2 c. 
in Figure 4. These models cover the full range of velocity pro- 
files and Lya emissivities. The model number associated with 
a given row shown in the top right corner of the right panels. 
The lower horizontal axis on the left column denotes the usual 
frequency variable x, which is converted to a velocity offset (in 
the Lya emitter's rest-frame) from the line center, on the upper 
axis. Note that positive values of x correspond to higher en- 
ergies, or bluer photons. The vertical axis on the left side has 
units which are such that the area under the curves corresponds 
to (27r)"' ( as in Fig. 1). The vertical axis on the right indicates 
the physical flux density in units of p Jy (ignoring scattering 
in the IGM, see § 8.1). The right panels panels show the sur- 
face brightness profiles. The horizontal axis denotes the impact 
parameter, measured in ". The units of the surface brightness 
on the vertical axis are ergs s"' cm"^ arcsec"^. The thick dot- 
ted horizontal line is the detection limit in the recent survey by 
Matsuda et al. (2004). 

First, we focus on the spectra which show that the Lya radiation 



escapes the collapsing protogalaxy predominantly blue shifted 
for all models. The extended models are not as blue as the cen- 
tral models, and they contain a faint red peak. We quantify the 
prominence of the red peak through the ratio [B] /[/?], defined 
as the number of photons in the blue peak over the number of 
photons in the red peak. 

Second, we focus on the surface brightness distribution. Model 
1. has a very flat profile, F{9) cx 0~", with a ^ 0.5. Model 2. has 
a steeper profile with a ^ 1 .0 and the central surface brightness 
is increased above Matsuda et al. (2004)'s detection limit. The 
profile that would emerge in the absence of scattering is over- 
plotted as the dashed curves. These curves illustrate the extent 
to which the surface brightness distribution becomes shallower 
because of scattering. The profiles for models 1 c. and 2 c. are 
more centrally concentrated {a ^ 2), which can be viewed sim- 
ply as an extension of the trend seen in the profiles of models 
1. and 2.: A steeper emissivity profile translates to a steeper 
surface brightness profile. 

6.2. Impact of Model Parameters on the Emitted Lya Features 

Here we discuss the impact of varying the model parameters 
described in § 5.2 on on the following properties of the emer- 
gent spectrum; (i) net blueshift of the blue peak, (ii) FWHM of 
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the blue peak, (iii) the location of the red peak (if any), (iv) the 
ratio and (v) the surface brightness profiles. 

6.2.1. Spectral Shape: Peak Morphology-The Blue Peak 
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Fig. 5. — This figure shows the morphology of the blue peak as afunction of 
various model parameters. The location of the maximum of the blue peak (tri- 
angles) and its FWHM (circles) as a function of 3 parameters for the extended 
(upper panel) and central (lower panels) models. The left panel, middle panel 
and right panel shows the dependence on a, Vamp and the optical depth r. The 
optical depth r is normalized to t\ , which is the optical depth from center to 
edge of model I. (see text). These trends suggest degeneracies will exist be- 
tween the quantities r, a and Vamp - for example, larger infall velocities can 
be mimicked by larger column densities and vice versa (and result in a larger 
blueshift of the Lya line). 



In Figure 5. we show the frequency offset of the blue peak's 
maximum (triangles) and its FWHM (circles) as a function of 
three parameters: a (left panels), Vamp ( middle panels) and the 
optical depth r from the center to the edge of the object (right 
panels). The upper and lower row of panels represents extended 
and central models, respectively. 

The left panels show that for the extended models, the de- 
pendence of the location of the maximum on a is weak, while 
its FWHM is slightly more affected by a (which is not surpris- 
ing, as we notice that the blue tail is more pronounced in the 
spectrum of model 2. than in model 1. in Fig. 4). On the other 
hand, the dependence on a is very strong for the central cases 
(which is also evident from the lower two rows in Figure 4). 
The central panels show that for the extended models the FWHM 
depends more strongly on Vamp than the peak's location. For the 
central models, the dependence of both quantities on Vamp is 
weaker than their dependence on a. Also, the dependence of 
the FWHM on Vamp for the central models is shghtly weaker 
relative to that of their extended counterparts. 
The right panels of Figure 5 shows the location of the maxi- 
mum of the blue peak and its FWHM for three different optical 
depths, T. The optical depth is normaUzed to that of model 1., 
denoted by ri. All models are identical in their dark matter 
and gas properties, but the Lya absorption cross section is de- 
creased and increased by a factor of 10 for the point with r/ri 



= 0.1 and 10.0, respectively. The figure demonstrates clearly 
that the blue peak's morphology is strongly affected by the 
value of T. Both the FWHM and velocity shift of the Lya line 
increase with r. These trends suggest that degeneracies will 
exist belwccn the values of t, Vamp and a - for example, larger 
infall velocities can be mimicked by larger column densities 
and vice versa. 

6.2.2. Spectral Shape: Peak Morphology-The Red Peak 




Fig. 6. — The location of the maximum of the red peak (-circles) and the 
minimum in the spectrum (-triangles) as a function of Vamp. The dotted line 
corresponds to the curve v = Vamp. Despite the fact that the included models 
vary in mass, redshift and a, the correlation is obvious, which suggests it is a 
generic property of the extended models. The location of the red peak poten- 
tially measitfes the gas infall speed. 

Figure 6 shows the location of the dip (circles) and the max- 
imum of the red peak (triangles) as a function of the amplitude 
of the velocity field for various models that vary in mass, red- 
shift and their assumed velocity field. The dotted line denotes 
V = Vamp. The exact parameters used for these models used in 
the plot can be read off from Table 2. 

The relation between Vamp and the location of the peak and 
dip is evident. Beyond Vamp ~ 360 km s"' the correlation breaks 
down. The reason is that the red and the blue peak merge for 
large Vamp (see Fig 5 and Fig 7 and their discussion below). The 
correlation holds for a wide range of models, which we take as 
evidence for it being a generic property of the extended models. 

As mentioned in § 6. 1 the blue to red ratio, [B] / [R] quantifies 
the prominence of the red peak. First we demonstrate how the 
transition occurs from a double peaked symmetric spectrum, 
known from the static cases, to the asynunetric spectra, that are 
associated with bulk motions. The left and right panel of Fig- 
ure 7 contains the (normalized) spectrum of models 1. and 1 
c, respectively, as the black-solid lines. The red-dotted lines 
correspond to the same models, in which the amphtude of the 
infall velocity has been reduced to 0.25vcirc. The blue dashed 
lines correspond to Vamp = (actually, Vamp 
for V 



amp 



10~^Vcirc, because 

: there is no emission for the extended cases accord- 
ing to eq. 11). 

The slight asynmietry that is present in both static cases, is 
intriguingly enough, due to deuterium. Deuterium also leaves a 
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sharp feature at jc = in the static, extended, case. The effect of 
deuterium on the other spectra was found to be negligible. Fig- 
ure 7 shows that increases strongly with Vamp for both 
the central and extended models. In the central cases, the blue 
and red peak remain well defined and separated as the ampli- 
tude of the infall velocity is increased. However, the amplitude 
of the red peak falls quickly to be undetectable. For the ex- 
tended cases, the FWHM of the blue peak increases, and that of 
the red peak decreases, until the peaks merge. The minimum in 
the spectrum that separates the blue and the red peak moves to 
redder frequencies with increasing Vamp- 

Figure 8 shows the value of as a function of Zvir 

(upper-right panel), v^^p (lower-right panel), a (lower-left panel) 
and Mtot (upper-left panel). The circles refer to the extended 
models, whereas the open stars refer to the central models. 
The upper-left panel shows that increases strongly with 

mass. The upper- right panel shows that decreases strongly 

with Zvir- In this panel only, we included the impact of the IGM, 
because the IGM introduces a strong z-dependence. A discus- 
sion of how the IGM was incorporated is delayed to § 8.1. The 
lower-right panel shows that the value of increases with 

the amplitude of the velocity field Va^p. For sufficiently large 
Vamp, the red peak is absorbed in the blue peak (see Fig 5), 
in which case — > oo. The lower-left panel shows that 
for the extended models is relatively insensitive to the 

value of a. For the central models, no photons emerged on the 
red side of the line for a < 0.25, causing oo. 
These results suggest that the secondary red peak may be de- 
tectable in the spectra of (proto) galaxies, provided the Lya 
emissivity is spatially extended. Furthermore, the red peak may 
be more easily detectable towards higher redshift. If detected, 
the location of the red peak gives a potentially useful measure 
of the gas infall speeds. 

6.2.3. The Surface Brightness Profile 

Finally, we focus on the impact of varying the input parame- 
ters on the emerging surface brightness profile. Figure 9 shows 
the surface brightness distribution as a function of Mtot, ct and 
j(r). The upper and lower panels are associated with extended 
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Fig. 7. — The effect of bulk motions in the gas on the emerging Ly« spec- 
trum is shown for the extended (left panel) and central (right panel) cases. 
The panels contain normalized spectra for three models that only differ in the 
amphtude of the velocity field. The black-solid line in the left/right panel cor- 
responds to model 1.11 c. respectively. The red dotted lines and blue dashed 
line correspond to the same models with Vamp reduced to 0.25vci,c and 75 c.) 
and Vamp = 0, respectively. The red peak becomes smaller with increasing infall 
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Fig. 8. — The ratio of the number of photons in the blue and red peak, 
as a function of various parameters. The fiducial model is model 
1. (see Table 2). Variation with: upper left panel: mass, upper right panel: 
redshift (In this panel only, we include the scattering in the IGM, since it intro- 
duces a strong redshift dependence, see § 8.1), lower right panel: t'amp, lower 
left panel: a. The red peak's prominence is enhanced towards higher redshifts, 
lower masses and lower infall velocities. 



and central models, respectively. Each panel contains three 
models for one particular mass, Mtot, which is given in the top 
right corner in units of 1O"M0. In all panels, the lowermost 
and uppermost curve corresponds to a = 1.0 and -0.5, respec- 
tively (also see Fig. 4), the middle curve has a = 0.25. Figure 9 
shows that the surface brightness profile flattens with increasing 
mass and a. 

Figure 10 shows the surface brightness profiles associated 
with model 1. (left panel) and model 2. (right panel). The 
two profiles are constructed from the reddest 15% (red solid 
line) and the bluest 15% (blue dotted line) of all Lya photons. 
The bluest 15% is more concentrated, especially for model 2., 
which has a = -0.5. This suggests that steepening of the sur- 
face brightness profile towards bluer wavelengths within the 
Lya line may be useful to diagnose gas infall. The amount 
of steepening may constrain the infall velocity profile. 

6.3. Summary of Intrinsic Lya Properties 
The results obtained above can best be summarized as: 

• Result 1. For all models, the Lya radiation emerges 
predominantly blueshifted, typically by 500-3000 km 
s"'.This can be explained as follows: 1) Photons prop- 
agating outwards in a collapsing gas cloud experience 
a net blueshift, because of energy transfer from the in- 
falling gas to the photons"* 2) Blueshifted Lya photons 
that are propagating outwards appear even bluer (and 
thus further in the wing) for the infalling gas, which fa- 

According to equation (3), Va ■ kj is negative (t^o and kj are anti-parallel for 
photons that are propagating outwards), while Va ko on average, which 
results in Jto > Xj. 
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Fig. 9. — The surface brightness as a function of angular distance in " from 
the center of the hya image. The upper and lower panels represent extended 
and central models, respectively. Each panel corresponds to a different mass, 
which is given in units of lO' 'Mq in the upper right comer of each panel. The 
velocity profiles have slopes a = —0.5 (steepest curve), 0.25 and 1.0 (shallow- 
est curve). The figure shows that the surface brightness profile becomes flatter 
with increasing mass. For a given mass, the exact slope is determined by how 
centrally concentrated the Lya emissivity is. 




8 
Distance (") 

Fig. 10. — This figure shows the frequency dependence of the surface bright- 
ness profile for the extended models. The surface brightness profile for the 
bluest {dotted blue line) and reddest {red solid line) 15% of the photons are 
shown. The left and right panel corresponds to models 1. and 2., respectively. 
The figure shows that for collapsing gas clouds, the bluest Lyo photons may 
emerge more centrally concentrated than the reddest Lya photons. 



cilitates their escape, while the reverse is true for red- 
shifted photons. These two effects yield a suppressed 
red and enhanced blue peak. 

Result 2. The spectra of the central models are system- 
atically bluer than those of the extended models. This 
result reflects that Lya photons that are inserted deeper 
inside the collapsing gas cloud, obtain a larger average 
blue shift as they propagate outwards. That is, they 



have to traverse a larger column of neutral infalling gas, 
which allows more energy to be transferred from gas to 
photons, which results in a larger average blue shift. 

• Result 3. The 'extended' models typically have an ad- 
ditional faint red peak in the spectrum separated by a 
minimum in the spectrum (the 'dip'). Both the loca- 
tion of the maximum of the red peak and the dip are 
uniquely determined by the infall velocity of the out- 
ermost gas layers. The reason for this is that some of 
the red Lya photons that are propagating outwards are 
exactly at resonance in the frame of the infalling gas, 
which significantly reduces their escape probability and 
produces the observed dip. The dip in the spectrum is 
therefore strongly linked to the gas' infall velocity (at 
the outermost gas shells). In § 8.2 we point out that de- 
spite the simplified nature of our models, in which the 
only bulk motions in the gas are radial and always in- 
flowing, this result is unlikely to change for more real- 
istic, more complex gas distributions. 

• Result 4. The Lya properties of the 'central' models 
depend more sensitively on a than those of the 'ex- 
tended' models, with the average blueshift of the spec- 
trum strongly increasing with decreasing a. This re- 
flects that all Lya photons in the central models are 
exposed to scattering by the gas in the central regions. 
Changing a from 1.0 to -0.5, enhances Vbuik by a fac- 
tor of (rvir/ r)'-^, which obviously becomes very large at 
small radii (r <C ryn). Due to the enhanced infall veloc- 
ity at small radii, more energy can be transferred from 
the gas to the photons in models with low values of a, 
which results in a larger average blue shift of the Lya 
line. 

• Result 5. The prominence of the red peak, quantified 
by the ratio [B] /[/?], increases with decreasing infall 
velocity, mass and total Lya optical depth. The ratio 

decreases with redshift for fixed other model 
parameters. These trends are easily understood: The 
larger the infall velocity, and the total Lya optical depth, 
the more energy is transferred from gas to photons as 
the photons work their way outwards. This results in a 
more prominent blue peak and larger [B] /[/?]. As ex- 
plained above, increasing a reduces the number of very 
blue photons and, subsequently, [B]/[^]. The increase 
of [B]/[^] with mass is caused by the increase of both 
the infall velocity, and the total Lya optical depth with 
mass. The decrease of [B]/[./?] with redshift caused by 
the increased probability with redshift for photons on 
the blue side of the line center to be scattered out of our 
line of sight in the IGM. This is discussed in more detail 
in§8.L 

• Result 6. The surface brightness profiles for the ex- 
tended models are typically flat, dF/dlnO = -1 to -0.5, 
with an increasing flatness with increasing mass and a. 
The same trends are seen for the central models, which 
consistently have steeper surface brightness distributions 
with dF/dln ^ -2. One reason that the surface bright- 
ness profiles of the extended models become flatter with 
increasing a is that the intrinsic emissivity becomes flat- 
ter with increasing a (see Fig 4). This is not the only 
reason though, since the same trend is seen among the 
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central models, which have identical j(r) (apart from the 
amplitude). As was mentioned above, for lower values 
of a the number of very blue Lya photons is enhanced, 
which increases the steepness of the surface brightness 
profile (see result 7. and its explanation below). The 
flattening of the surface brightness profile with mass is 
caused by: 1) The projected emissivity on the sky be- 
comes increasingly shallow with mass, which is due to 
the mass dependence of the virial radius ^ (for the ex- 
tended models). 2) The Lya optical depth increases with 
mass, which reduces the probability that photons prop- 
agate through the outer gas layers without being scat- 
tered. 

• Result 7. A flattening of the surface brightness profile 
towards longer wavelengths within the Lya line is in- 
dicative of infall. This can be explained as follows: the 
bluest photons emerge mainly from the center (as ar- 
gued above) and are unlikely to get scattered in the outer 
gas layers. This is because these photons appear even 
bluer in the frame of the infalling gas, and thus further 
in the line wing. The reverse is true for the reddest pho- 
tons, which appear much closer to the line center in the 
frame of the infalling gas and therefore likely scattered 
in the outermost gas layers. This produces a more dif- 
fuse appearance of the reddest photons compared to that 
of the bluest. 

7. FLUORESCENT LYa EMISSION FROM SELF-IONIZED 
PROTOGALAXIES 

In this section, we focus on the observable properties of flu- 
orescent Lya, such as its surface brightness distribution and its 
spectrum. For this purpose, an ionizing point source was in- 
serted at the center of our sphere. Its total flux in ionizing pho- 
tons is Lion (in ergs s*') and its spectrum is a power law of the 
form Jiy^r) = Ji){r){i^ / i^iY^ . The slope (3 depends on the na- 
ture of the ionizing sources. We choose /3 = 1.7, for the point 
source to represent a quasar. The value of Joir) is obtained 
by Jo(r) = (/?- l)Lion/[16fL7r^r^]. The standard equations of 
photoionization equilibrium were solved numerically with the 
code. The optical depth to the first shell was assumed to be 0. 
This was used to calculate its neutral fraction of hydrogen. This 
allows the optical depth to the next shell to be computed, which 
also determines its neutral fraction. This procedure is repeated 
until the outer shell is reached. The resulting changes in the 
intrinsic Lya spectrum are presented below. 

7.1. Impact of a Central Ionizing Source on the Properties of 
the Emerging Lya Emission 

In the presence of an ionizing source we set the Lya emis- 
sivity to be given by: 

j(r)=J^aaB[l-XH(r)Yn^if(r) hvo+fc ^ — (14) 

where is the average number of Lya photons produced per 
recombination, which is ^ 68% for case B and ^ 42% for case 
A recombinations. Case A/B applies when the gas is optically 
thin/thick to all Lyman series photons, respectively. The sec- 
ond term is the same as used for the extended cooling models 

'if the projected emissivity were plotted as a function of O/Oy^^y in which ^yir 
is the angle the virial radius subtends on the sky, it would be identical for all 
masses. 



(eq. 11). The main reason for keeping this term is to make 
the transition from the previous extended models to the current 
models smooth. 

The gas in model 7. becomes fully ionized at Lion = 10''^ 
ergs s"'. We caution that this estimate ignores small-scale gas 
clumping. If we denote the clumping factor by C = (p^) / {p)^, 
then this transition to full ionization occurs at ^ C times higher 
ionizing luminosity (which would allow for C times higher 
Lya luminosity) 

Increasing Lion beyond this value, merely reduces the residual 
neutral fraction, and consequently, A'hi- The 'phase transition' 
from neutral to fully ionized dramatically affects the emerging 
Lya spectrum. In Figure 1 1 we show the Lya spectrum and 
surface brightness distribution for four values of Lion: Lion = 
{black-solid line), Lion = 5 x 10^^ ergs s~' (red-dotted line). 
Lion = 10"*^ ergs s~' (blue-dashed line) and Lion = 10"*^ ergs s~' 
(green-long-dashed line). The results are presented differently 
than above: The left panel shows the surface brightness distri- 
bution rotated by 90°. The right panels show spectra of these 
models at four different impact parameters, 9, indicated by the 
horizontal dotted lines. For example, the lowest panel shows 
the spectrum emerging from 6 ~ 0-2", the second lowest panel 
shows the spectrum emerging from 6^2- 4", etc. The main 
reason for presenting the data in this fashion is to study the vari- 
ations in the spectra across the object. Note that the spectra are 
arbitrarily scaled for visualization purposes. 

We note that in the the absence of an ionizing source, the 
spectrum (black-solid line) is fairly constant across the object, 
the central region being sUghtly bluer. For low ionizing lumi- 
nosities. Lion = 5 X lO''^ ergs s"\ photo-ionization enhances the 
Lya emission from the central regions. Enhanced emission in 
the center enhances the number of photons in the bluest parts of 
the spectrum (§ 6.3), and steepens the surface brightness pro- 
file (observation 6. in § 6.3). These are indeed the trends that 
can be seen in Figure 11. When Lion is increased by a further 
factor of two, the gas sphere as a whole is fully ionized, which 
produces several effects. The overall blueshift of the spectrum 
(blue-dashed line ) almost completely vanishes, and its FWHM 
is significantly smaller. The main reason is the reduced column 
of neutral hydrogen, which reduces the importance of broaden- 
ing of the Lya line by resonant scattering. Instead, the broad- 
ening is determined by the bulk velocity of the gas, despite the 
fact that each Lya photon scatters lO'' times before it es- 
capes. The reason that the Lya spectrum is barely affected by 
scattering, is that each Lya photon mainly scatters near the site 
of its creation, where the photon is at resonance. A chance en- 
counter with a rapidly moving atom can put the photon in the 
wing of the Une, which allows it to escape in a single flight 
from the local resonant region (in which the optical depth is 
To ^ lO-'-well within the regime of moderate opacities, Adams 
1972). These photons escape this region in a typical symmetric 
double peaked profile around the local Une center. Velocity gra- 
dients in the gas facilitate the photon's escape from the rest of 
the object, effectively without 'seeing' shells at other radii. The 
frequency with which the photon finally emerges at the virial 
radius therefore does contain information on the bulk velocity 

*We mentioned in § 1 that the total Lya emissivity can be greatly boosted 
by embedded ionizing sources. A requirement for this boost is that the gas is 
clumpy, which increases the recombination rate and allows a larger fraction of 
the ionizing radiation to be reprocessed into Lya photons, rather than escaping, 
e.g. (Partridge & Peebles 1967; Haiman & Rees 2001 ; Alam & Miralda-Escudfi 
2002). Because dumpiness is not included in our model, this boost in the 
Lya luminosity is not achieved. 
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Fig. 1 1. — This figure displays the Lya surface brightness distribution {left panel) and spectra emerging {right panels) from a collapsing gas cloud containing 
an ionizing source for four different ionizing luminosities. The spectra in the different panels are centered on the impact parameter indicated by the black-dotted 
lines. Model 1. is represented by the black— solid line (Ljo,, = 0). The red-dotted line, blue— dashed line and green-long-dashed line corresponds to model 1. plus 
an additional central ionizing source with Li„„ = 5 X lO"*- ergs s"' , L,„„ = 10"*^ and lO"** ergs s"' , respectively. As Ljo^ increases, the spectral line becomes narrower 
and more symmetric around the line center. 



of the site where the photon was generated. The double peak in 
the emission line in regions 1-111 is a consequence of this local 
resonant scattering. The main reason for the enhanced promi- 
nence of the dip at large impact parameters is that the velocity 
gradient projected along the line of sight is smaller there, which 
increases the total column of hydrogen a Lya photon effectively 
sees. 

The reason that the spectrum from region IV is narrower is 
that this spectrum mainly contains photons that emerge from 
the central regions, in which the bulk infall velocity is low (a = 
1.0). The effect of the velocity field and Lion ion the emerging 
spectrum is investigated in more detail in § 7.2. 

7.2. Constraints from Lya Emission on Gas Kinematics in the 
Presence of a QSO 

In the previous section we argued that the gas kinematics can 
be constrained by the Lya emission despite the possibility that 
the Lya photon scatters up to ^ 10"* times. In this section we 
show that, although the amplitude of the bulk velocities in the 
gas can be constrained well by the gas, the sign is much harder 
to determine, since the emerging spectra for infall and outflow 
are almost identical for fully ionized gas. In this section, we 
ignore contribution from cooling to the emissivity; the second 
term on the RHS in equation (14) is set to 0. The reason is that 
we want to explore the similarities of the effects of scattering. 



so we prefer to adopt the same emissivity profiles. 

In Figure 12 we show spectra at four different impact param- 
eters for three different models. The blue-dashed lines repre- 
sent model 7. containing a central ionizing source with Lion = 
lO"*"* ergs s"' (as in Fig. 1 1). These spectra are compared to the 
no-scattering profiles, which refer to the profiles that would be 
observed if the gas were transparent to Lya (shown as black- 
solid lines). The red-dotted lines represent a model in which 
the sign of Vamp is reversed. Each panel contains the spectra at 
the impact parameter indicated by the roman numerals in the 
top-right corner (as in Fig. 11). The only prominent differ- 
ence between all spectra, is that the models in which scattering 
occurs are shifted by several Doppler widths relative to the no- 
scattering spectrum. Infall and outflow models are slightly blue 
and redshifted shifted, respectively. The reasons for this shift 
have been discussed in § 6.3. 

Although not shown in this paper, we found the surface bright- 
ness distributions for the optically thin infall and outflow mod- 
els to be the identical (also see ZM02, who found the same). 
This leads to the conclusion that for gas that is only moder- 
ately thick to Lya scattering, tq <^ lO"*, the Lya spectrum can 
be used to determine its velocity dispersion, despite the fact that 
Lya photon up may scatter up to 10^- 10"* times. However, the 
sign of the velocity field can only be determined when the ex- 
act Lya line center is known, to an accuracy of ^ 10 km s"' 
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Fig. 12. — This figure shows the Lyo spectra emerging from optically thin 
outflowing and infalling gas. Each panel corresponds to the impact parameter 
indicated by the roman letter (as in Fig. 11). The spectra shown by the blue- 
dashed lines represent a photo-ionized infall model, while red-dotted lines 
represent the outflow models in which the sign of 

^'amp is reversed. The black- 
solid lines represent both infall and outflow models in the absence of scattering. 
The figure shows that when the gas is optically thin to Lya photons, tq < 10^, 
infall and outfiow models can produce nearly identical Lya spectra, except for 
a small shift in frequency. 



(which is probably not possible in practice). This introduces an 
ambiguity in the interpretation of the Lya blob found around a 
z = 3 quasar by Weidinger et al. (2004, 2005). Although they 
attribute the observed spectra and surface brightness profiles to 
gas infall, our results suggest that outflowing gas can produce 
similar characteristics. We will explore this issue in more depth 
in paper 11. 

8. DISCUSSION 
8.1. The IGM Transmission 

As mentioned in § 5.2, we have focussed on the intrinsic Lya 
emission from collapsing protogalactic gas clouds. The subse- 
quent Lya transfer through the IGM, which is defined to lie 
beyond the virial radius, has been ignored so far Even though 
the IGM is almost fully ionized at most redshifts considered in 
this paper, a small residual neutral fraction is enough to scatter 
some fraction of the Lya photons from our line of sight. For 
an IGM that is comoving with the Hubble flow, the mean trans- 
mission is 1 and (e~^) redward and blueward of the line center, 
respectively (e.g. Madau 1995). Note that (e""^) is to a good ap- 
proximation constant, apart from the small-scale fluctuations in 
wavelength caused by the Lya forest (i.e., overall redshift evo- 
lution is not significant). The exact value of (e""^) is a function 
of redshift and is well known from studies of the Lya forest 
(e.g. Songaila & Cowie 2002; Lidz et al. 2002; McDonald et al. 
2001) and is tabulated in Table 1. 

The naive way to incorporate the impact of the IGM, is to 
give a photon that escaped with x>Q and jc < a weight (e""^) 



Table 1 
Transmission of IGM 



2.4 



3.0 



3.8 



4.5^ 5.2- 



5.7- 



6.0^ 



< > 0.8 0.68 0.49 0.25 



0.09 
0.02 



0.07 
0.003 



0.004 
0.003 



(1) From McDonald et al. (2001) 

(2) From Songaila & Cowie (2002) 



and 1, respectively. This is how the dependence of on 
redshift was calculated in § 6.2.2 (see Fig. 8). 

In § 6 we expressed the total Lya luminosity, LLya, as a func- 
tion of our model parameters. To correctly convert this to a flux 
on earth, denoted by fhya, is not straightforward, due to the 
frequency dependence of the IGM transmission. Motivated by 
the observation that the Lya emerges mainly blueshifted from 
our models, we estimate the total detectable Lya flux from: 
fhya ~< e"'^ > Li^ya/4-TTdl. We obtain: 



/i 



Lya 



10- 



lOiiM, 







1.75 



1.2 



x(t^)'''^z<4.0 



ergs 



(15) 



where we used < e"^ >= 0.68[(1 +z)/4]-'^ for z < 4 and < 
>= 0.031[(l+z)/7]-** ' for z > 4, where the last approxi- 
mation was taken from (Songaila & Cowie 2002). Note that 
this approximation fails at z > 6, where the mean transmission 
drops even more rapidly ^. We caution that the strong break 
in the slope of the redshift dependence of the flux (eq. 15) at 
z > 4 should not be taken literally. In reality, the variation of 
the this slope with redshift should be smoother. The main point 
of eq. (15) however, is that because the IGM becomes increas- 
ingly opaque towards higher redshifts at z > 4, cooling radia- 
tion of a fixed luminosity is more difficult to detect at z ^ 4 
than at z ^ 4. Quantitatively, the total detectable Lya cool- 
ing flux does not scale with d^^ (which can be approximated as 
oc (1 +z)~^'^^), but as cx< e""^ > dj^ (which can be approximated 
as c>c (1 +z)~'') for z > 4.0. We note that in the the proximity of 
bright QSO's, significantly more flux may be transmitted (e.g. 
Cen & Haiman 2000; Madau & Rees 2000; Cen et al. 2005), 

The above model of IGM transmission is appropriate for an 
IGM that is comoving with the Hubble flow. The gravitational 
pull exerted by the massive objects considered in this paper, 
causes their surrounding gas to decouple from the Hubble flow 
and move inwards. As a consequence, photons that escape our 
object with a small redshift, can Doppler shift into resonance 
in this infalling IGM, which increases their probability of being 
scattered out of the line of sight. The result of this is that the 
IGM may erase a fraction of the red part Lya fine (Barkana & 
Loeb 2004; Santos 2004). We postpone a detailed discussion of 
this effect to Paper 11, but note that it barely affects the quoted 
values of [B] in Figure 8. 

8.2. Robustness of our Results 

'To allow ftya to be written as a simple function of (1 + z), we approximated 
the luminosity distance by dtij,) = 3.6 X 10''([l+z]/5)"^'* Mpc, which is within 
< 6% of its actual value, in the range z = 3 — 10. 
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The simplified nature of our models clearly introduces un- 
certainties to the results presented in this paper. Nevertheless, 
several results obtained in this paper should be robust. The con- 
clusion that Lya emerges blue-shifted from a collapsing cloud 
will not change for different gas distributions, since it simply re- 
hes on Lya photons encountering gas with a bulk velocity that 
opposes their propagation direction. The only requirement is 
that the Lya photons do meet gas on their way out. If the gas is 
very clumpy, and only covers part of the sky from the site where 
the Lya is generated, then it could escape without a net blue 
shift. For the massive protogalactic halos considered in this pa- 
per the covering factor is expected to be well above one, and we 
do not expect Lya to escape without encountering any gas. As 
argued by Birnboim & Dekel (2003), virial shocks may be ab- 
sent in halos with Merit ^ 2 x 10' ^Mq, which suggests that even 
if the gas is clumpy, in the protogalactic phase, the clumps are 
not embedded in a hot tenuous medium for masses below Merit- 
Although the mass of our fiducial model lies slightly above this 
quoted critical mass, its exact value depends on redshift, metal- 
Ucity of the gas, etc. In this case the gas is also optically thick 
to Lya photons between the denser clumps, making our calcu- 
lations more appropriate. 

The presence of the dip in the emergent spectrum redward of 
the hne center (Result 3. in § 6.3) is most hkely not an artifact 
of the simplified nature of our models. The representation of 
the gas dynamics in our models, as a smooth and radial infall, 
is a clear oversimplification. In reality, gas infall is likely to be 
clumpy, and may have significant tangential velocity compo- 
nents. However, even for this type of complex gas distribution 
and velocity field, a the dip is Ukely to be present, since it only 
requires Lya photons redward of the line center to be at reso- 
nance in the frame of the gas in the outermost shells. Since the 
location of the dip will depend on the local velocity of the gas in 
the outermost shells, the position of the dip would vary across 
the object for more complex gas motions. However, as long as 
the gas cloud as a whole is collapsing, the overall escape prob- 
ability for Lya photons redward of the line center is reduced, 
which may show up as a noticable (but broader) dip. Observa- 
tional support for this claim is that Barkana & Loeb (2003) have 
found evidence for the presence of a dip in the Lya spectra of 
several quasars, which is believed to be caused by infalling gas. 

Similarly, that the bluest Lyaphotons can freely stream out- 
wards through a colunan of collapsing gas, while the reddest 
Lya photons get scattered in the outermost regions, should be 
valid regardless of the exact gas distribution. It is not clear 
however, what fraction of centrally produced Lya escapes. The 
presence of even a small amount of dust is sufficient to ab- 
sorb Lya photons (Neufeld 1990; Chariot & Fall 1991). If the 
medium is clumpy, then Lya photons may still escape by scat- 
tering of (dusty) cold clouds embedded in a tenuous medium 
(Neufeld 1991; Haiman & Spaans 1999; Hansen & Oh 2006). 
The red peak, which is produced further out, suffers less from 
dust absorption. The net effect of destroying preferentially the 
bluer photons would weaken the spectral dependence of the sur- 
face brightness distribution for the extended models. There- 
fore, although a fiattening of the surface brightness distribution 
with increasing wavelength (within the Lya line) is indicative 
of gas infall, its absence, however, does not argue against in- 
fall. Once detected, the amount of flattening could constrain 
the infall velocity profile. Detailed results, such as the exact 
shapes of the spectra, are very likely to change for reaUstic gas 
configurations. For reUable detailed predictions, we should ap- 



ply the code to the output of cosmological simulations which 
can provide the Lya code with realistic density, velocity and 
temperature fields (as in Tasitsiomi 2005). However, the cal- 
culations presented in this paper relating to peak morphology 
do offer simple insights in what may be plausible Lya proper- 
ties of collapsing protogalaxies. In Paperll, we demonstrate the 
usefulness of our models, by applying them to several observa- 
tions. 

The gas temperature in this paper was assumed to be lO'* 
°K. Turbulent motions in the gas may raise the effective tem- 
perature of the gas, which was the reason ZM02 used T ~ 
10^ °K. The exact temperature of the gas affects the frequency 
redistribution of Lya photons undergoing resonant scattering 
and therefore affects the results in which the broadening and 
blueshift of the Lya emission line are dominated by resonant 
scattering. However, the temperature dependence of our results 
is weak. This is best illustrated for the case of a static, uniform 
sphere, in which the emergent spectrum is entirely determined 
by resonant scattering. If all Lya photons are injected in the 
center of the sphere at x = 0, the emergent spectrum peaks at 
Xp = ±0.92(aro)'/^ (§ 2). For a fixed total column of hydrogen 
from the center to the edge of the sphere, the velocity offset 
from the line center of the two peaks depends on temperature 
as: Vshift oc J'/*. A change in temperature by a factor of 10 
shifts our lines by a factor of ;^ 1 .5. 

8.3. Comparison with Previous Work 

In this paper we used a Monte Carlo technique that is similar 
to those used in previous papers (e.g. Ahn et al. 2000, 2002; 
Zheng & Miralda-Escude 2002; Cantalupo et al. 2005). The 
following discussion briefly summarizes how our work differs 
from and complements the work previously done. Ahn et al. 
(2002) and Ahn et al. (2003) have focussed primarly on static, 
extremely thick slabs which represented star burst galaxies and 
optically thick expanding superwinds. The Lya from their su- 
perwind models emerged with a net redshift, for the same rea- 
son our models predict blueshifted Lya emission. Zheng & 
Miralda-Escude (2002) concentrated on fluorescent Lya emis- 
sion from Damped Lyman a (DLAs) systems illuminated by 
an externally generated UV background (either diffuse, or pro- 
duced by a nearby quasar). They found self-shielding to pro- 
duce a core in the surface brightness profiles, which differs 
from the surface brightness profiles of any of our models. Ad- 
ditionally, the spectra emerging from their fluorescent DLAs 
are double peaked and symmetric around the line center. How- 
ever, this symmetry is mainly due to the fact that Zheng & 
Miralda-Escude (2002) choose their gas cloud to be static. Can- 
talupo et al. (2005) also considered fluorescent Lya radiation 
from high column density systems using high resolution hydro- 
simulations. They found the symmetry in the spectra of fluores- 
cent Lya radiation to vanish for more reaUstic velocity fields, 
e.g. collapsing gas clouds produce a more prominent blue peak, 
similar to the spectra presented in this paper. However, fluore- 
sent emission in response to the diffuse ionizing background 
only reaches flux levels of 10~^^ ergs s~' cm~^, which is well 
below most of the signals calculated in this paper. The amount 
of fluorescent Lya emission can be boosted significantly in the 
presence of a nearby bright quasar, which should be easily iden- 
tifiable. It appears possible to observationally distinguish be- 
tween our models and those discussed above. 

9. SUMMARY AND CONCLUSIONS 
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We presented Lya Monte Carlo radiative transfer calcula- 
tions through neutral (§ 5), and self-ionized (§ 7), collapsing, 
spherically symmetric gas clouds. To provide a test case for 
our code that is tailored to spherical symmetry, we derived a 
new analytic solution for the Lya spectrum emerging from a 
uniform, static sphere (§ 3.2, Appendix C). We found that a 
cosmological abundance of deuterium regulates the transfer of 
Lya photons over a narrow frequency range (corresponding to 
40-50 km s-i at T = lO'* °K, see § 4). Although the results in 
the present paper are not affected by the presence of deuterium, 
we point out it may affect the analysis for systems with hydro- 
gen column densities in the range 10^^ — 10^" cm"-^, for T = 10* 
°K. 

We studied the effects of gas kinematics and Lya emissivity 
as a function of radius, j(r), on the Lya spectrum emerging 
from the collapsing clouds and on the Lya surface brightness 
distribution, for a range of models with different masses and 
redshifts (§ 6.2), which are summarized in Table 2 

The emergent Lya spectrum is typically double-peaked and 
asymmetric. In practice, we find that the blue peak is signifi- 
cantly enhanced, which results in an effective blueshift of the 
Lya line. This blueshift may easily be as large as ~ 2000 km 
s"' (§ 6.2.1). The total blueshift of the line increases with the 
total column density of neutral hydrogen of the gas and with its 
infall speed (§ 6.2.1, Fig 5). The prominence of the red peak 
is enhanced towards higher redshift, lower mass and lower in- 
fall speeds (6.2.2, Fig 8). If detected, the shift of the red peak 
relative to the line center may give a potentially useful mea- 
sure of the gas infall speed (§ 6.2.2, Fig 6). Furthermore, a 
steepening of the surface brightness profile towards bluer wave- 
lengths within the Lya line is a useful diagnostic for gas in- 
fall. The amount of steepening may constrain the infall veloc- 
ity profile (§ 6.2.3, Fig 10). A more complete summary of all 
results involving the surface brightness profiles, the morphol- 
ogy of the blue and red peak, and their relative prominence is 
given-together with a qualitatitive explanation-in § 6.3. 

The 'fluorescent' (or 'self-ionized') cases (§ 7) relate to mod- 
els of collapsing gas clouds with an embedded ionizing source. 
As its ionizing luminosity increases, the Lya line becomes nar- 
rower and more symmetric around the Une center (Fig 11). When 
the gas is fully ionized, we found the line widths in the spec- 
trum to be a good measure for the velocity distribution of the 
gas, despite the fact that in some cases each Lya photon may be 
scattered up to 10^ - 10"* times. We caution that for these opti- 
cally thin cases, infall models have almost the same spectrum as 
outflow models (and identical surface brightness profiles, also 
see ZM02). The only difference is that for infall/outflow the 
emission Une is shifted in frequency to the blue/red, respec- 
tively (Fig. 12), by an amount which decreases with the bright- 
ness of the central ionizing source. We will explore this issue 
in more depth in paper II. 

The effective blueshift of the Lya emission from collapsing 
protogalaxies has profound consequences for its detectability, 
as radiation blueward of the Lya line is subject to IGM ab- 
sorption (§ 8.1), an effect that increases rapidly with redshift 
{e~'^ cx (1 +z)~** for z > 4). We expect the Lya fluxes from high 
redshift (z > 4.0) collapsing protogalaxies to be very weak, un- 
less the intrinsic Lya luminosity is very high (eq. 15). In Paper 
II, where we use the main results of this paper to interprete ob- 
servations of Lya emitters, the impact of the IGM on our results 
is discussed in more detail. 

Ongoing and future surveys of Lya emitters are providing us 



with useful constraints on the early stages of galaxy formation. 
We hope the work presented in this paper can contribute to the 
interpretation of recent and of forthcoming observations of Lya 
emitting galaxies. 
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APPENDIX A 

THE MODELS 

The table below contains all models that we studied in this 
paper. For each model, its associated set of model parame- 
ters Mtot, Zvir, 01 and Vamp is given. The infall velocity, Vamp has 
units of the circular velocity for (given in the fourth column). 
The first column contains the model number. Two models that 
only differ in their choice of j{r) have the same model num- 
ber, but the 'central' (see the discussion above) model has an 
additional label 'c' attached . For reference, the last colunm 
contains the figure numbers in which results pertinent to each 
model are shown. 
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APPENDIX B 
TESTS OF THE CODE 



B.l. Testing Core Scattering 

The test discussed in § 3.2 mainly addresses scattering of photons in the wing, since > 99% and > 99.98% of the photons that 
emerge from the optically thick sphere are in the wings, for tq = 10^ and 10^, respectively. To calculate these numbers, the transition 
from core to wing was defined at |x| = 3. Therefore, the emerging spectrum is insensitive to the performance of the code for scatterings 
in the core. Note that this is exactly why it is allowed to deploy the accelerated scheme (§ 3.1.1) to calculate emerging Lya spectra 
from optically thick media. To test the accuracy of the code for scatterings in the core we compare the redistribution function, 
qn(x,x') (§ 2), as extracted from the code, with the theoretical solutions given by Hummer (1962) and Lee (1974) (their qnb(x,x')). 
In the upper right panel of Figure 1 we show the cases with x = 0,x = 2 and x = 5. The agreement between the histograms generated 
by the code, and the theoretical curves is excellent. 

An additional test for core scatterings involves the well studied problem of the mean number of scatterings a Lya photon under- 
goes, N, before it escapes from the slab (see Adams 1972). To calculate the average number of scatterings a Lya photon undergoes 
in our code, we use the weighting scheme described by Avery & House (1968). This weighing procedure reduces the average number 
of scatterings by a factor of ^ 2. The main reason is that the distribution of the number of scatterings has a long tail to large number 
of scatterings, which increases N. The weighting scheme described by Avery & House (1968), gives less weight to the photons that 
are in the high-AT tail, which brings the number of scatterings of individual photons closer to the mean. The dots in the lower left 
panel denote as a function of tq as extracted from our simulation, whereas the solid line is the theoretical prediction by Harrington 
(\912>)N= 1.61to. The agreement at large tq is impressive, but breaks down when tq ^ 10^, which is in the transition regime between 
moderate and extreme optical depths. 



B.2. Testing the Code when the Gas is not Static 

To test the vaUdity of the code in a case when the gas has bulk motions, we compare our results with LR99, who used a Monte Carlo 
method to calculate the Lya transfer through a fully neutral, homogeneous, uniform intergalactic medium that undergoes Hubble 
expansion. They demonstrated that ionizing sources prior to reionization, are embedded in large Lya halos that are redshifted relative 
to the host galaxy. Initially, resonant scattering moves the Lya photons into both the red and blue wing of the line profile, which 
dramatically increases their mean free path (as discussed in § 2). The increased mean free path enables Hubble expansion to redshift 
the Lya photons during their flight. The blue photons return to the line core, while the red photons become even redder. The Lya 
photons that have redshifted far enough from the line center stream freely to the observer. LR99 use a simple Monte Carlo approach 
to calculate the spectrum (their Fig. 2) and surface brightness profiles (their Fig. 1) of these Lya halos. We quantify this below: 

A photon initially at frequency x (in this case x <C 0), propagating through an expanding IGM, redshifts in proportion to its traversed 
distance, dx = -H(z)ds/vth- Here, H(z) is the Hubble parameter at redshift z, ds is the differential traversed physical distance, and 
Vth is the thermal velocity of the atoms in the IGM. As the photon redshifts, its mean free path increases, and consequently Hubble 
expansion redshifts the photon even further, etc. The maximum optical depth a photon initially at frequency x can see through a 
uniform IGM is given by: 



Tr = 



: / ds nH(s)(7(s) = hhcto f ds(p(x [s]) = -^=^ / ds—r— = / —r= (Bl) 

For T = 10°/r,z=10, and the WMAP cosmological parameters. 



_ -3.8x103 _ X. 

Tx « = — , (B2) 

X X 

in which is the 'critical' frequency, which marks the frequency beyond which the neutral IGM is optically thin. Photons with 
X <C X* (note that x is negative ) stream freely through the neutral IGM towards the observer. LR99 define their frequency variable in 
terms of this critical frequency as: x = x/x* . Note that blue- and redshifted photons have a negative and positive x, respectively. 

We caution that the results presented in LR99 are not exactly adequate to test our code in the case of bulk motions. In LR99 the 
Lya photons redshift during their flight, while the energy before and after scattering is equal (i.e. Xg = x,). In the case the shells 
have bulk motion, the photon's energy changes in a scattering event according to (eq. 3), while the photon's energy remains constant 
during its flight. For example: A photon is emitted at x = and r = in the direction k. Consider the situation in which each scattering 
leaves ^unchanged. In LR99, the photon continuously redshifts in between fully coherent scatterings. In the situation in which the 
gas shells are radially expanding (according to Hubble's law), the photons frequency will remain the same, x = 0. Note that the 
mean free path of this photon is the same in both situations. This photon emerges at a hypothetical boundary in the same number of 
scatterings in both cases, but in LR99 it is redder. The situation sketched above is very rare, but it is easy to convince oneself that 
any photon, along an arbitrary trajectory, emerges sUghdy redder in LR99's scenario. In practice, this difference is small, but leaves 
a noticeable difference in the spectrum. 

LR99 clearly describe their Monte Carlo simulation, which can easily be reproduced and adequately modified to describe Lya 
transfer through gas radially expanding outwards from r = according to the Hubble law, Vbuik =H r, where H is the Hubble constant. 
The calculated spectrum for this modified LR99 simulation is slightly bluer than their original spectrum shown in their Figure 2. This 
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modified curve is shown in the lower right panel of our Figure 1 as the green-solid line. Note that LR99 used dimensionless units, 
with all cosmological parameters scaled out. The histogram shown was computed with our own code. Lya photons were injected in 
the center of a sphere in which the shells are expanding outwards according to Hubble' s law, Vbuik = H{z) r. In our code, the variables 
are not dimensionless. We pick Zs =10 (which is the physical model LR99 translate their results to) and the density of the gas shells 
is uniform and «// ^ 2 x 10"^ (1 +Zsf' cm"-'. The gas temperature was set to 10 K (a = 1.5 x 10"^, although the exact temperature 
barely matters). To facilitate comparison with LR99n, this plot uses their frequency variable, i ~ -3.8 x 10^ x. As the figure shows, 
the agreement is excellent. 



APPENDIX C 

ANALYTIC SOLUTION FOR A STATIC UNIFORM SPHERE 

In this Appendix we derive an analytic solution for the Lya spectrum emerging from a uniform, static sphere of gas, similar to 
the well known solution for a plane-parallel slab derived by Harrington (1973) and Neufeld (1990). Following Harrington (1973); 
Neufeld (1990), our starting point is the radiative transfer equation after having made the Eddington approximation. In spherical 
coordinates this reads: 

d^J(r,x) 2 dJ(r,x) ^ ,2, x 2,,, x o n 2., n 
~d?~ r~dT~ " (x)Ko(Jir,x)-3(j)(x)K^ir), 

L{r)= [ (j){x!)q{x,x!)J{x! ,r)dx! + (i){x) — . (CI) 
i-00 47r 

Rybicki & dell' Antonio (1994) showed that the integral operator can be written as a differential operator, after applying the Fokker- 
Planck approximation: 



f 

J— a 



' 1 d dJ(x r) 

(f){x)q{x, x)J{x , r)dx ^ (/>(x)7(x, r) + -— (j)(x) — . (C2) 

2 ox ax 



After defining the variable cr, by dx/da = (3 /2)^/^(j>{x) (note that in the line wing, (j)(x) = a/[-K^l^y?\ and therefore a = (27r/27)i/ V/a), 
the transfer equation simplifies to: 

d^J{r,x) 2 dJ(r,x) . d^J{r,x) _ ^,2.. 27W ™. 

Following Urmo (1955), this equation is solved using the Sturm-Liouville expansion which consists of an eigenfunction expansion of 
the form: 



Jir,x) = J2En(r)h„ix). (C4) 



Let each eigenfunction, E„(r), obey the homogeneous differential equation: 

d'^E„ir) 2dE„{r) ^ sinAr cosAr 

a 2 + 5 — + A^£„(r) = E„{r)=C— — +D— — , (C5) 

or^ r or \r \r 

where C and D are arbitrary constants. Note that to prevent the solution from diverging at r = 0, D = 0. The value of A is obtained 
from the boundary condition. The flux can be written as (e.g. Rybicki & Lightman 1979; Unno 1955; Harrington 1973; Neufeld 
1990, although the first authors have an extra factor of tt in the definition of the flux F): 

F{r,x) = - — . (C6) 

iKo<p{x) or 

At the boundary, where the radius is denoted hy r = R, all photons are propagating outwards, i.e. no photons get scattered back 
through the surface. In this case, the flux F(R,x) = 2J(R,x). For this equaUty b(r) = was assumed. Each function E„{r) should obey 
the boundary conditions, which gives the following constraint on A: 

AcotA/g = -^(/.(x) + ^^A„ = ^7r«( 1- ^ +0[l/i3<t>ix)KoR+2f]] , n=l,2,3,... (C7) 

2 R " \ 3<p{x)K.i)R + 2 I 

The sequence of solutions for A was obtained after making the same approximations as in Unno (1955); Harrington (1973). Solutions 
with different values of n are orthogonal. Each £'„(r) can be normalized by multiplying the function by A„/-\/27r/?. Then Jir,x) is: 

rr ^ i. ^ ^ ^i"'^"'" ^ sin(A„r) 

J(r,x) = y C„ h„{x) — = ^== > ^ h„{u). (C8) 

, A„r \/2-kR , r 



20 



We substitute this back into the radiative transfer equation (eq. C3). Following Harrington (1973), the function 2>(jp-{x) (which is 
peaked strongly around x = 0) is replaced by \/65(cr). The solution for any particular h„{a) can be be found by multiplying the 
radiative transfer equation on both sides by the m* eigenfunction, sin(Amr)/r\/27r/? followed by integration over J A-Kp-dr. 



Q, 



-f 

Jo 



Anrdr 



2 sin(A„0 . 



Kr). 



The solution for h„{(j) is given by: 



, , . V6ko „ exp(-A„|a|/Ko) 



which in combination with equation (C8) completes the solution for J{r,x). 

C. 1 . Analytic Solution for a Central Source. 
First consider the situation in which all photons are emitted at radius r^, then: 

]{r) = 6{r-rs), Qm = , 

when substituted into equation (C8) we get: 



sin(A„r) / \/6ko Kcp(-A^^|aj/Ko)\ 

Sim ^ j 



n=l * 



V6 



16tt^R rr. 



y^sin(AnOsin(Anr,) 



n=l 



exp(-An|g|/Ko) 

An 



The spectrum at the boundary is obtained by setting r = R. From equation (C7) we know: 



sin(A„/J) = cos(A„7?) 



\nR 



(-1)"; 



An 



Therefore, 



Using 



we find: 



J{R,a): 



1-1 .5(j)(x)RKo 1 .5(/>(x)ko 

•\/6 ■^2\/7rx^^ /— nTT |c7| \ sinA„rj 



167r2/?2 



> -5 (-l)"exp — 

3aKo \ R Kny r, 

n=\ 



n=l 



V(-1)V=-; , for|x|<l ,cosx 

2V7rjc2 / 1 



smx = 



2i 



coshx : 



e^ + e 



JiR,a) = 



32iTT^R\ 3aKo Vi+exp(i|^ + j7rr,//?) l+exp(^-j7rr,//?)' 
\/6 / sin(7rrs//?) 
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^ u 

l+exp(^-j7rr,//?)^ 



■)■ 



(C9) 



(CIO) 



(Cll) 



(C12) 



(C13) 



(C14) 



(C15) 



(C16) 



'«««^os(W^) + cosh(^^S)^ 
After r, and multiphcation by 47r/?^ (to obtain the total emerging flux density at the surface) we obtain the equation given in 
equation (9) (setting Rkq = tq): 



Jix) = 



lAaTQ Vi+cosh 



l22lL\i\ 
27 aro 



Solution for sphere, 



which can be compared to the previously known analytic solution for a slab (Harrington 1973; Neufeld 1990): 



V6 



24V^«^oVi+cosh 



54 OTo 



Solution for slab. 



(C17) 



(C18) 
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APPENDIX D 

DEVIATIONS FROM PARTIAL COHERENCE 

In § 2 we mentioned that for all applications discussed in this paper, the scattering is partially coherent. The main reason for this 
is that the Einstein A coefficient is very high, which translates to a very short spontaneous decay time, fspon = = 10"** s. Because 
the hydrogen atoms spend such a short time in the n = 2 state, the chance the atom gets perturbed, while in this excited state, is tiny. 
Below this is quantified: 

D. 1 . Collisions with Electrons 

Collisions between the excited hydrogen atom and electrons, can either cause the energy of the re-emitted Lya photon to differ 
from that of the incoming photon, or result in the loss of the Lya photon. The first process occurs when a collision puts the atom in 
any state from which it radiatively decays back into the 2p state. The last process occurs when the collision puts the atom in a state 
from which the 2p state is not accessible via a radiative transition. An example of such a state is the 2s level. Osterbrock (1962) 
gives the collision rate from the 2p to 2s state, at 7" = lO"' °K , C2p2s ~ 2 x W^ng s~\ where is the number density of electrons. 
This collision rate dominates the collision rates into all other states by several orders of magnitude (e.g Giovanardi et al. 1987; Chang 
et al. 1991). Therefore, only the collision rate from 2p to 2s is relevant for our purposes. Given that the maximum densities in our 
models reach rie ^ 10^ cm"^, the rate at which excited hydrogen atoms get perturbed is at most 10"' s~\ which is nine orders of 
magnitude below the spontaneous decay rate, and therefore negligible. 

D.2. Absorption of Additional Photons 
The Lya scattering rate for a given hydrogen atom is, 

Pa=A-K I dv ; S « — ; i^o ~ 5 X 10 JalX S , (Dl) 

J hv hvQ c 

where we assumed that J(x) = J a is constant and wrote J a = Ja,i\ x 10"^' ergs s"' cm~^ sr"' Hz"'. 7^21 is highest in the central 
regions in cases in which all Lya is inserted at r = 0. An estimate of 7„,2i can be obtained from: 

where /S.v is the line width of the Lya line, which we set conservatively to IQ^'^Ui). Our lines are typically much broader, which 
would reduce This implies a scattering rate of ~ 10"^ s"' at r = 1 kpc. Since the Einstein B-coefficient for stimulated emission 
equals that of absorption, the rate at which a hydrogen atom in the excited state absorbs Lya photons equals that of a hydrogen atom 
in the ground state. The rate at which stimulated emission of Lya occurs per hydrogen atom is therefore > 12 orders of magnitude 
lower than A21 for r > 1 kpc, and can be safely ignored. 

In the presence of a continuum source , the photo-ionization rate per hydrogen atom from the 2p state by any photon with £ > 3.4 
eV, emitted by a central source with a powerlaw spectrum Jiu) oc (for hv>2)A eV) is 

Pion = 47r/ dv ; s ' wSx 10 — ^ rr Lo,o-7^ ^ ' (^3) 

J34ev V 10^^ ergs s ' / \/5 + 2.75/ V r I 

where an^iv) is the photo-ionization cross section from the « = 2 state at frequency v, which is ~ 1.4 x 10"'^(i^/fHa)~^'^^> where 
hiAia = 3.4 eV (e.g. Cox 2000, p. 108). The total luminosity of the source is given by Ltot- Photo-ionization out of the n = 2 state by 
continuum photons emitted by the central source is clearly negUgible. 

The final mechanism by which the atom may be perturbed in its excited state is by absorbing other line photons. Since the Einstein 
A coefficients are weaker for all other transitions in hydrogen, the photo-excitation cross section to other levels is less than the Lya 
absorption cross section. For the same number density of these other line photons, the photoexcitation rate is less than the Lya 
scattering rate calculated above, which was negligible. Therefore, also this process can be safely ignored, which completes our 
justification for only considering partially coherent scattering (as compared to incoherent scattering). 



